Percolation, Morphogenesis, and Burgers Dynamics in Blood Vessels Formation 
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Experiments of in vitro formation of blood vessels show that cells randomly spread on a gel matrix 
autonomously organize to form a connected vascular network. We propose a simple model which 
reproduces many features of the biological system. We show that both the model and the real 
system exhibit a fractal behavior at small scales, due to the process of migration and dynamical 
aggregation, followed at large scale by a random percolation behavior due to the coalescence of 
aggregates. The results are in good agreement with the analysis performed on the experimental 
data. 
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The problem of morphogenesis and self-organization in 
biological systems is both an actual and a long standing 
one In living beings, complex structures arise from 

the aggregation of separate constituents. Such struc- 
tures, although irregular, often show some sort of hidden 
order, like self-similarity and scaling laws [HQ. Scal- 
ing laws are known to emerge not only in the physics 
of phase transition [f| and in percolation 0, but also 
from several kinds of aggregation dynamics, and are of- 
ten considered as fingerprints of the process which led 
to the formation of the structure itself 0, 0- A typ- 
ical example of natural structure characterized by non 
trivial scaling laws is the vascular network In re- 
cent years many experimental investigations have been 
performed on the mechanism of blood vessel formation 
101. Cells are cultured on a gel matrix and their migra- 
tion and aggregation observed through videomicroscopy. 
Tracking of individual trajectories shows marked persis- 
tence in the direction, with a small random component 
superimposed [Tlj . The motion is directed towards zone 
of higher concentration of cells, suggesting that chemo- 
tactic factors play a role. Cells migrate over distances 
which are an order of magnitude larger than their radius 
and aggregate when they get in touch with one of their 
neighbors. In a time of the order of 10 hours they form a 
continuous multicellular network which can be described 
as a collection of nodes connected by chords. The mean 
chord length is approximately independent on the initial 
cell density n, with an average value I = 200 ± 20 ^m for 
n between 100 to 200 cells/mm 2 . Observing the formed 
structures, one finds that, by varying the initial cell den- 
sity, there is a percolative transition. Below a critical 
value n c ~ 100 cells/mm 2 groups of disconnected struc- 
tures are formed (Fig.QJi). For n > n c a single connected 
network is instead visible (Fig. |IJ>,c). For higher values 
of n one observes a sort of "swiss cheese" pattern (Fig. 
Eli). In the biological system the percolating property is 
of physiological relevance, since it is directly linked with 



the functionality of blood vessels. 

In this paper we propose a theoretical model which 
turns to be in good agreement with experimental obser- 
vations. It allows to well reproduce both the observed 
percolative transition and the typical scale of observed 
vascular networks. With respect to standard percolative 
models, here an important role is played by migration 
and dynamical aggregation of particles. The model ap- 
pears thus as a possibly new and physically interesting 
representative in the class of percolation models. We 
characterize the phase transition from the point of view 
of scaling laws, both in the model and in the real system, 
computing critical indices and the fractal dimension of 
the percolating cluster through extensive numerical simu- 
lations and analysis of the experimental data. The model 
describes the cell population as a continuous density field 
n and velocity v; it also assumes the presence of a con- 
centration field c of soluble factor. Cells are modeled as a 
fluid accelerated by gradients of the soluble factor. The 
latter is supposed to be released by cells, diffuse, and 
degrade in finite time, in agreement with experimental 
observations. These assumptions give rise to the follow- 
ing equations: 
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where D, a, r, are respectively the diffusion coefficient, 
the rate of release and the characteristic degradation time 
of soluble mediators, and fi measures the strength of cell 
response. Initial conditions are given in the form of a set 
of randomly distributed bell-shaped bumps in the den- 
sity field, representing single cells initially spread in the 
system, with zero velocities and zero concentration of 
the soluble factor. In Eqs. Q a multidimensional Ric- 
mann Ijlbjl equation is coupled to the diffusion equation 



pcjl. while (|la(l is just a continuity equation expressing 
conservation of the total cell number. Multivaluedness 
of solutions to (|ib|) can be prevented by adding a vis- 
cous term V 2 v in the r.h.s. One than gets the analog 
of Burgers' equation, which has been used to describe 
growth of molecular interfaces 01 an d the emergence of 
network-like patterns in the large-scale distribution 
of masses of the Universe 0, 0|. Instead of viscos- 
ity, we introduce a phenomenological, density dependent 
pressure term — V0(n), in the r.h.s. of (|lb|l . with <p(n) 
zero for low densities, and rapidly increasing for densi- 
ties above a threshold er -2 , where a = 30 /jm is the mean 
cell radius, in order to describe the fact that cells do not 
interpenetrate. 

The model provides a natural lenght scale ro = V Dt, 
which is the effective range of the interaction mediated by 
the soluble factor. This suggests that, starting from the 
above mentioned initial conditions, Eqs. Q should de- 
velop network patterns characterized by the same scale. 
Direct measurements of D and r give D ~ 10~ 7 cm 2 /s, 
t ~ lh, from which one gets ro ~ 200 /im, a value which 
is in good agreement with the mean chord length t ex- 
perimentally observed. 

We performed numerical simulations of model on 
square boxes of size L = 1, 2, 4, 8 mm, with periodic 
boundary conditions, using a finite volume method |l6| . 
We used the experimental values for D and t, while 
the unknown parameters /i, a were set to one (which 
amounts to an appropriate rescaling of time and the con- 
centration field c). The initial conditions were given by 
throwing the same number of cells as in the biological ex- 
periments in random positions inside the box, with zero 
velocities and zero concentration of the soluble factor, 
with a single cell given initially by a Gaussian bump of 
width a and unitary weight in the integrated cell density 
field n. Starting from these condition, Eqs. were nu- 
merically integrated. The simulation was stopped when 
the vascular network was formed, or a stationary state 
was reached. We observed a remarkable correspondence 
between simulations and experiments both in the the dy- 
namic evolution details (not shown) and in the reproduc- 
tion of the percolating patterns (Fig. |2J . 

In order to study the percolative transition we com- 
puted between 100 to 200 realizations of the system for 
each box size and cell density, with different initial po- 
sitions of the cells. The simulation gives as a result a 
continuous density field n. In order to study percolation 
properties we partitioned the box in little square boxes of 
size 2 -6 mm, forming the sites of a square lattice, and set 
each site as filled if the mean cell density inside the cor- 
responding little box was larger than the threshold cj~ 2 , 
empty otherwise. We then identified clusters of nearest 
neighbor filled sites. In Fig. the results of four sim- 
ulations, for different values of the density, are shown. 
Black areas represent regions filled with cells, that is re- 
gions where density exceeds the threshold, white areas 




FIG. 1: Experimental pictures of vascular networks, on a sys- 
tem of size L = 2 mm, obtained starting from four different 
values of the initial cell density: a) 50 cells/mm 2 ; b) 100 
cells/mm 2 ; c) 200 cells/mm 2 ; d) 400 cells/mm 2 . 




FIG. 2: Numerical simulations, on a system of size L = 2 mm, 
obtained starting from four different values of the initial cell 
density: a) 50 cells/mm 2 ; b) 100 cells/mm 2 ; c) 200 cells/mm 2 ; 
d) 400 cells/mm 2 . 

represent the underlying substrate. As in the case of ex- 
perimental data, a percolative transition is observed at 
some critical cell density n c . Furthermore the model re- 
produces quite well the typical structure of the vascular 
network, with chords of length t ~ 200 \xm for a wide 
range of cell densities, in good agreement with the ex- 
pected value ro- 

For each set of realizations, corresponding to a value 
of the box size L and of the mean cell density n, we 
identified clusters of nearest neighbor filled sites, and 
measured the following quantities: 1) the percolation 
probability II, defined as the fraction of realizations in 
which a percolating cluster appears; 2) the infinite clus- 
ter density P, defined as the density of sites belonging 
to the infinite (greatest) cluster; 3) the mean cluster size 
S = (1/N) J2 S n sS 2 , where ./V is the total number of sites 
and n s is the number of clusters of size s, excluding from 
the sum the infinite cluster. 
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FIG. 3: Data collapse of a) percolation probability II and b) 
mean cluster size S, for the numerical simulations with system 
sizes L = 1, 2, 4, and 8 mm. 



In percolation models, in presence of a second order 
transition at some mean cell density n c , these quantities 
show a characteristic singular behavior near the perco- 
lation point. Namely, the percolation probability is zero 
below the transition and one above it; the infinite cluster 
density, which plays the role of order parameter, vanishes 
as \n — n c \P from above; and the mean cluster size has 
a singularity \n — n c | -7 on both sides of the transition. 
For finite size systems, the same quantities obey near the 
transition to the following relations: 
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Using Eq. (|2a() . one can find the critical density n c as 
the point where the curves II(n, L) for different sizes L 
cross. The estimated value is n c — 94 ± 1 cells/mm 2 . 
Then, the critical index v is found by plotting H(n, L) 
as a function of (n — n c )L 1 / l/ , where n c is held fixed to 
the value previously computed, and looking for the value 
of v that gives the best data collapse. The error on the 
computed value of the index is estimated looking for the 
range of values that gives an acceptable data collapse. 
Holding the value of n c and of v fixed, the other two 
relations Eq. (|2b|) and (|2cl) are used to estimate the value 
of the indices (3 and 7, by looking for the data collapse 
of l// v P(n, L) and L^l" S(n, L) as a function of (n - 
7i c )L 1 / l/ . We observed that while n c is sensitive to the 
choice of the threshold in a neighborhood of the natural 
value cr -2 , the critical exponents do not depend on it. 
In Fig. 03 the data collapses for II and S are shown. The 
estimated values of the critical indices are reported in the 
second column of Tab. B] 

Within the errors, the values obtained are the same 
that one expects for the simple random percolation 
model, that is when the sites of the lattice are uncorre- 
cted, and each one occupied with the same probability. 
As critical indices are linked to the large scale structure of 
the percolating cluster, this means that, on such a large 
scale, the structure of the vascular network is mainly de- 
termined by the initial random positioning of the cells, 



critical index 


model 


experiments 


random percolation 


V 


1.33 ±0.08 


not measured 


1.333 


7/1/ 


1.83 ±0.05 


1.78 ±0.12 


1.792 


P/v 


0.11 ±0.01 


not measured 


0.104 


D 


1.87 ±0.03 


1.85 ±0.10 


1.896 



TABLE I: Critical indices measured on the numerical model 
and on experimental data, and compared with the exact val- 
ues of random percolation. 



and is not altered by the dynamical process of migration 
and aggregation. On the other hand, on smaller distances 
clusters are quite different from the ones of random perco- 
lation. For example the mean fraction of occupied sites, 
which is a quantity sensible to the small scale structure 
of clusters, is about 0.2 at the percolation point, much 
lower than the corresponding figure for random perco- 
lation, that is 0.59. A quantity that can give us infor- 
mations about the structure of the percolating cluster at 
different scales is the density p(r) as a function of the 
radius. This is defined as the mean density of sites be- 
longing to the percolating cluster, inclosed in a circle of 
radius r centered at one site belonging to the cluster, and 
averaged over different centering sites and different real- 
izations of the system. For a fractal object with fractal 
dimension D, this should scale as p(r) ~ r D ~ d . For the 
percolating cluster of random percolation at the critical 
point, one expects a fractal dimension D = 1.896. We 
computed p(r) for our model at the percolation point 
n = n c (see Fig. 0J, and found two different regimes de- 
pending on the scale. For r > r c , with r c = 0.77 ± 0.08 
mm, a behavior compatible with random percolation is 
found, with fractal dimension D = 1.87 ± 0.03. On the 
other hand, for r < r c a different exponent is found, 
D = 1.50 ± 0.02. This exponent may be the signature 
of the dynamic growth process driven for r < r c by the 
rapidly oscillating components of the concentration field. 

We have repeated the measures described above on a 
set of experimental data, consisting of 28 digital pho- 
tographs (with a resolution of 1024 2 pixels) of mature 
structures obtained on a matrix of size L — 2 mm, 
starting from initial densities of 50, 75, 100, . . ., 200 
cells / mm 2 . The photographs are the results of two exper- 
imental sessions, each performed in duplicate, so that 4 
photographs are available for each density value. The 
data show the presence of a percolative transition at 
n c ~ 125 cells/mm 2 , and are compatible with the values 
7/1/ = 1.83 obtained from numerical simulations. The 
available data are not enough to obtain v and (3/u with 
a reasonable precision. Estimated values of the critical 
indices are reported in the third column of Tab. [I] Fi- 
nally, we measured the density of the percolating cluster 
as a function of the radius at the percolation point, and 
found a behavior surprisingly similar to the one found in 
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FIG. 4: Density of the percolating cluster as a function of the 
radius for numerical simulations and experimental data. 



the numerical model (see Fig. 0J: for r > r c the fractal 
dimension of the percolating cluster is D = 1.85 ± 0.10 
which is compatible with D = 1.87±0.03 from the numer- 
ical simulations and with the value D = 1.896 of random 
percolation. For r < r c instead D = 1.48 ± 0.05. 

In conclusion, we have introduced a model which de- 
scribes the formation of the vascular network with a typ- 
ical scale lenght ro, and is in good agreement with exper- 
imental data. The percolative phase transition appears 
to be second order, at least to the accuracy allowed by 
the limitations in size which are intrinsic to our numeri- 
cal experiments. Comparison of critical indices suggests 
that the phase transition falls in the universality class of 
random percolation, even in presence of migration and 
dynamical aggregation. This is confirmed by the frac- 
tal dimension of the percolating cluster, on scales larger 
than r c ~ 0.8 mm, that is in agreement with the one 
expected for random percolation. On the other hand we 
found that, if observed on scales smaller than r c , the vas- 
cular network shows a different scaling behavior, with a 
fractal dimension D ~ 1.5, that may be the signature of 
the dynamical process that led to the formation of the 
network. 

Thus, the model appears to be quite successful in de- 
scribing in vitro experiments, where all the parameters 
are under control and one can easily tune the cell den- 
sity. It would also be interesting to understand the rele- 
vance of this picture to vascular networks in living beings. 
It is known [lOj that during embriogenesis blood vessels 
form through local (in situ) differentiation of a fraction 
of endothelial cells, that afterwards assemble in a vascu- 
lar labyrinth, a process similar to the one described in 
this paper. However, in this case it is more difficult to 
have complete control over all the relevant parameters. 



In particular, one would expect that the cell density is 
fixed through some feedback mechanism, in order to max- 
imize efficiency Indeed, some experiments on 
in vivo capillaries show compact space-filling structures, 
with fractal dimension D = 2, which would indicate that 
cell density is tuned by Nature to a value larger than 
the percolation threshold. More experiments, focusing 
on the relevance of the percolation mechanism, and the 
nature of cell density self-regulation, are needed to shed 
light on these issues. 
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